home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / clang / pgp20src.zip / MPILIB.C < prev    next >
C/C++ Source or Header  |  1992-08-05  |  51KB  |  1,474 lines
  1. /*    C source code for multiprecision arithmetic library routines.
  2.     Implemented Nov 86 by Philip Zimmermann
  3.     Last revised 27 Nov 91 by PRZ
  4.  
  5.     Boulder Software Engineering
  6.     3021 Eleventh Street
  7.     Boulder, CO 80304
  8.     (303) 541-0140
  9.  
  10.     (c) Copyright 1986-92 by Philip Zimmermann.  All rights reserved.
  11.     The author assumes no liability for damages resulting from the use 
  12.     of this software, even if the damage results from defects in this 
  13.     software.  No warranty is expressed or implied.  The use of this 
  14.     cryptographic software for developing weapon systems is expressly 
  15.     forbidden.
  16.  
  17.     These routines implement all of the multiprecision arithmetic 
  18.     necessary for number-theoretic cryptographic algorithms such as 
  19.     ElGamal, Diffie-Hellman, Rabin, or factoring studies for large 
  20.     composite numbers, as well as Rivest-Shamir-Adleman (RSA) public 
  21.     key cryptography.
  22.  
  23.     Although originally developed in Microsoft C for the IBM PC, this code 
  24.     contains few machine dependancies.  It assumes 2's complement 
  25.     arithmetic.  It can be adapted to 8-bit, 16-bit, or 32-bit machines, 
  26.     lowbyte-highbyte order or highbyte-lowbyte order.  This version
  27.     has been converted to ANSI C.
  28.  
  29.  
  30.     The internal representation for these extended precision integer
  31.     "registers" is an array of "units".  A unit is a machine word, which
  32.     is either an 8-bit byte, a 16-bit unsigned integer, or a 32-bit
  33.     unsigned integer, depending on the machine's word size.  For example,
  34.     an IBM PC or AT uses a unit size of 16 bits.  To perform arithmetic
  35.     on these huge precision integers, we pass pointers to these unit
  36.     arrays to various subroutines.  A pointer to an array of units is of
  37.     type unitptr.  This is a pointer to a huge integer "register".
  38.  
  39.     When calling a subroutine, we always pass a pointer to the BEGINNING
  40.     of the array of units, regardless of the byte order of the machine.
  41.     On a lowbyte-first machine, such as the Intel 80x86, this unitptr
  42.     points to the LEAST significant unit, and the array of units increases
  43.     significance to the right.  On a highbyte-first machine, such as the
  44.     Motorola 680x0, this unitptr points to the MOST significant unit, and
  45.     the array of units decreases significance to the right.
  46.  
  47.     Modified 8 Apr 92 - HAJK
  48.     Implement new VAX/VMS primitive support.
  49.  
  50. */
  51.  
  52. /* #define COUNTMULTS */ /* count modmults for performance studies */
  53.  
  54. #ifdef DEBUG
  55. #ifdef MSDOS
  56. #include <conio.h>
  57. #define poll_for_break() {while (kbhit()) getch();}
  58. #endif    /* MSDOS */
  59. #endif    /* DEBUG */
  60.  
  61. #ifndef poll_for_break
  62. #define poll_for_break()  /* stub */
  63. #endif
  64.  
  65. #include "mpilib.h"
  66.  
  67. short global_precision = 0; /* units of precision for all routines */
  68. /*    global_precision is the unit precision last set by set_precision.
  69.     Initially, set_precision() should be called to define global_precision
  70.     before using any of these other multiprecision library routines.
  71.     i.e.:   set_precision(MAX_UNIT_PRECISION);
  72. */
  73.  
  74. #ifdef PORTABLE
  75. /*************** multiprecision library primitives ****************/
  76. /*    The following portable C primitives should be recoded in assembly.
  77.     The equivalent assembly primitives are externally defined under
  78.     different names, unless PORTABLE is defined.  See the header file
  79.     "mpilib.h" for further details.
  80.  
  81.     The carry bit mechanism we use for these portable primitives
  82.     will only work if the size of unit is smaller than the size of
  83.     a long integer.  In most cases, this means UNITSIZE must be
  84.     less than 32 bits -- if you use the C portable primitives.
  85. */
  86.  
  87. typedef unsigned long int ulint;
  88. #define carrybit ((ulint) 1 << UNITSIZE)
  89. /* ...assumes sizeof(unit) < sizeof(unsigned long) */
  90.  
  91. boolean mp_addc
  92.     (register unitptr r1,register unitptr r2,register boolean carry)
  93.     /* multiprecision add with carry r2 to r1, result in r1 */
  94.     /* carry is incoming carry flag-- value should be 0 or 1 */
  95. {    register ulint x;    /* won't work if sizeof(unit)==sizeof(long) */
  96.     short precision;    /* number of units to add */
  97.     precision = global_precision;
  98.     make_lsbptr(r1,precision);
  99.     make_lsbptr(r2,precision);
  100.     while (precision--)    
  101.     {    x = (ulint) *r1 + (ulint) *post_higherunit(r2) + (ulint) carry;
  102.         *post_higherunit(r1) = x;
  103.         carry = ((x & carrybit) != 0L);
  104.     }
  105.     return(carry);        /* return the final carry flag bit */
  106. }    /* mp_addc */
  107.  
  108.  
  109. boolean mp_subb
  110.     (register unitptr r1,register unitptr r2,register boolean borrow)
  111.     /* multiprecision subtract with borrow, r2 from r1, result in r1 */
  112.     /* borrow is incoming borrow flag-- value should be 0 or 1 */
  113. {    register ulint x;    /* won't work if sizeof(unit)==sizeof(long) */
  114.     short precision;    /* number of units to subtract */
  115.     precision = global_precision;
  116.     make_lsbptr(r1,precision);
  117.     make_lsbptr(r2,precision);
  118.     while (precision--)    
  119.     {    x = (ulint) *r1 - (ulint) *post_higherunit(r2) - (ulint) borrow;
  120.         *post_higherunit(r1) = x;
  121.         borrow = ((x & carrybit) != 0L);
  122.     }
  123.     return(borrow);        /* return the final carry/borrow flag bit */
  124. }    /* mp_subb */
  125.  
  126. #undef carrybit
  127.  
  128.  
  129. boolean mp_rotate_left(register unitptr r1,register boolean carry)
  130.     /* multiprecision rotate left 1 bit with carry, result in r1. */
  131.     /* carry is incoming carry flag-- value should be 0 or 1 */
  132. {    register short precision;    /* number of units to rotate */
  133.     register boolean nextcarry;
  134.     precision = global_precision;
  135.     make_lsbptr(r1,precision);
  136.     while (precision--)    
  137.     {    
  138.         nextcarry = (((signedunit) *r1) < 0);
  139.         /* nextcarry = ((*r1 & uppermostbit) != 0); */
  140.         *r1 <<= 1 ;
  141.         if (carry) *r1 |= 1;
  142.         carry = nextcarry;
  143.         pre_higherunit(r1);
  144.     }
  145.     return(nextcarry);    /* return the final carry flag bit */
  146. }    /* mp_rotate_left */
  147.  
  148. /************** end of primitives that should be in assembly *************/
  149. #endif    /* PORTABLE */
  150.  
  151.  
  152. /*     The mp_shift_right_bits function is not called in any time-critical
  153.     situations in public-key cryptographic functions, so it doesn't 
  154.     need to be coded in assembly language.
  155. */
  156. void mp_shift_right_bits(register unitptr r1,register short bits)
  157.     /*    multiprecision shift right bits, result in r1. 
  158.         bits is how many bits to shift, must be < UNITSIZE.
  159.     */
  160. {    register short precision;    /* number of units to shift */
  161.     register unit carry,nextcarry,bitmask;
  162.     register short unbits;
  163.     if (bits==0) return;    /* shift zero bits is a no-op */
  164.     carry = 0;
  165.     bitmask = power_of_2(bits)-1;
  166.     unbits = UNITSIZE-bits;        /* shift bits must be < UNITSIZE */
  167.     precision = global_precision;
  168.     make_msbptr(r1,precision);
  169.     while (precision--)    
  170.     {    nextcarry = *r1 & bitmask;
  171.         *r1 >>= bits;
  172.         *r1 |= carry << unbits;
  173.         carry = nextcarry;
  174.         pre_lowerunit(r1);
  175.     }
  176. }    /* mp_shift_right_bits */
  177.  
  178.  
  179. #ifndef VMS
  180.  
  181. short mp_compare(register unitptr r1,register unitptr r2)
  182. /*    Compares multiprecision integers *r1, *r2, and returns:    
  183.     -1 iff *r1 < *r2
  184.      0 iff *r1 == *r2
  185.     +1 iff *r1 > *r2
  186. */
  187. {    register short precision;    /* number of units to compare */
  188.  
  189.     precision = global_precision;
  190.     make_msbptr(r1,precision);
  191.     make_msbptr(r2,precision);
  192.     do    
  193.     {    if (*r1 < *r2) 
  194.             return(-1);
  195.         if (*post_lowerunit(r1) > *post_lowerunit(r2)) 
  196.             return(1);
  197.     } while (--precision);    
  198.     return(0);    /*  *r1 == *r2  */
  199. }    /* mp_compare */
  200.  
  201. #endif /* VMS */
  202.  
  203. boolean mp_inc(register unitptr r)
  204.     /* Increment multiprecision integer r. */
  205. {    register short precision;
  206.     precision = global_precision;
  207.     make_lsbptr(r,precision);
  208.     do    
  209.     {    if ( ++(*r) ) return(0);    /* no carry */
  210.         post_higherunit(r);    
  211.     } while (--precision);
  212.     return(1);        /* return carry set */
  213. }    /* mp_inc */
  214.  
  215.  
  216. boolean mp_dec(register unitptr r)
  217.     /* Decrement multiprecision integer r. */
  218. {    register short precision;
  219.     precision = global_precision;
  220.     make_lsbptr(r,precision);
  221.     do    
  222.     {    if ( (signedunit) (--(*r)) != (signedunit) -1 )
  223.             return(0);    /* no borrow */
  224.         post_higherunit(r);    
  225.     } while (--precision);
  226.     return(1);        /* return borrow set */
  227. }    /* mp_dec */
  228.  
  229.  
  230. void mp_neg(register unitptr r)
  231.     /* Compute 2's complement, the arithmetic negative, of r */
  232. {    register short precision;    /* number of units to negate */
  233.     precision = global_precision;
  234.     mp_dec(r);    /* 2's complement is 1's complement plus 1 */
  235.     do    /* now do 1's complement */
  236.     {    *r = ~(*r);
  237.         r++;    
  238.     } while (--precision);
  239. }    /* mp_neg */
  240.  
  241. #ifndef VAXC /* Replaced by a macro under VAX C */
  242.  
  243. void mp_move(register unitptr dst,register unitptr src)
  244. {    register short precision;    /* number of units to move */
  245.     precision = global_precision;
  246.     do { *dst++ = *src++; } while (--precision);
  247. }    /* mp_move */
  248.  
  249. #endif /* VAXC */
  250. void mp_init(register unitptr r, word16 value)
  251.     /* Init multiprecision register r with short value. */
  252. {    /* Note that mp_init doesn't extend sign bit for >32767 */
  253.     register short precision;    /* number of units to init */
  254.     precision = global_precision;
  255.     make_lsbptr(r,precision);
  256.     *post_higherunit(r) = value; precision--;
  257. #ifdef UNIT8
  258.     *post_higherunit(r) = value >> UNITSIZE; precision--;
  259. #endif
  260. #ifdef VAXC
  261.  
  262.     unitfill0( r, precision); /* Use character insts. on a VAX */
  263.  
  264. #else /* VAXC */
  265.  
  266.     while (precision--)
  267.         *post_higherunit(r) = 0;
  268.  
  269. #endif /* VAXC */
  270.  
  271. }    /* mp_init */
  272.  
  273.  
  274. short significance(register unitptr r)
  275.     /* Returns number of significant units in r */
  276. {    register short precision;
  277.     precision = global_precision;
  278.     make_msbptr(r,precision);
  279.     do    
  280.     {    if (*post_lowerunit(r)) 
  281.             return(precision);
  282.     } while (--precision);
  283.     return(precision);
  284. }    /* significance */
  285.  
  286.  
  287. #ifndef VAXC    /* Replaced by a macro under VAX C */
  288.  
  289. void unitfill0(unitptr r,word16 unitcount)
  290.     /* Zero-fill the unit buffer r. */
  291. {    while (unitcount--) *r++ = 0;
  292. }    /* unitfill0 */
  293.  
  294. #endif /* VAXC */
  295.  
  296. /* The macro normalize(r,precision) "normalizes" a multiprecision integer 
  297.    by adjusting r and precision to new values.  For Motorola-style processors 
  298.    (MSB-first), r is a pointer to the MSB of the register, and must
  299.    be adjusted to point to the first nonzero unit.  For Intel/VAX-style
  300.    (LSB-first) processors, r is a  pointer to the LSB of the register,
  301.    and must be left unchanged.  The precision counter is always adjusted,
  302.    regardless of processor type.  In the case of precision = 0,
  303.    r becomes undefined.
  304. */
  305.  
  306.  
  307. /* The macro rescale(r,current_precision,new_precision) rescales
  308.    a multiprecision integer by adjusting r and its precision to new values.  
  309.    It can be used to reverse the effects of the normalize 
  310.    routine given above.  See the comments in normalize concerning 
  311.    Intel vs. Motorola LSB/MSB conventions.
  312.    WARNING:  You can only safely call rescale on registers that 
  313.    you have previously normalized with the above normalize routine, 
  314.    or are known to be big enough for the new precision.  You may
  315.    specify a new precision that is smaller than the current precision.
  316.    You must be careful not to specify a new_precision value that is 
  317.    too big, or which adjusts the r pointer out of range.
  318. */
  319.  
  320.  
  321. /*    The "bit sniffer" macros all begin sniffing at the MSB.
  322.     They are used internally by all the various multiply, divide, 
  323.     modulo, exponentiation, and square root functions.
  324. */
  325.  
  326.  
  327. int mp_udiv(register unitptr remainder,register unitptr quotient,
  328.     register unitptr dividend,register unitptr divisor)
  329.     /* Unsigned divide, treats both operands as positive. */
  330. {    int bits;
  331.     short dprec;
  332.     register unit bitmask;
  333.     if (testeq(divisor,0))
  334.         return(-1);    /* zero divisor means divide error */
  335.     mp_init0(remainder);
  336.     mp_init0(quotient);
  337.     /* normalize and compute number of bits in dividend first */
  338.     init_bitsniffer(dividend,bitmask,dprec,bits);
  339.     /* rescale quotient to same precision (dprec) as dividend */
  340.     rescale(quotient,global_precision,dprec);
  341.     make_msbptr(quotient,dprec); 
  342.  
  343.     while (bits--)
  344.     {    mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
  345.         if (mp_compare(remainder,divisor) >= 0)
  346.         {    mp_sub(remainder,divisor);
  347.             stuff_bit(quotient,bitmask);
  348.         }
  349.         bump_2bitsniffers(dividend,quotient,bitmask); 
  350.     }
  351.     return(0);
  352. } /* mp_udiv */
  353.  
  354.  
  355. #define RECIPMARGIN 0    /* extra margin bits used by mp_recip() */
  356.  
  357. int mp_recip(register unitptr quotient,register unitptr divisor)
  358.     /* Compute reciprocal (quotient) as 1/divisor.  Used by faster modmult. */
  359. {    int bits;
  360.     short qprec;
  361.     register unit bitmask;
  362.     unit remainder[MAX_UNIT_PRECISION];
  363.     if (testeq(divisor,0))
  364.         return(-1);    /* zero divisor means divide error */
  365.     mp_init0(remainder);
  366.     mp_init0(quotient);
  367.  
  368.     /* normalize and compute number of bits in quotient first */
  369.     bits = countbits(divisor) + RECIPMARGIN;
  370.     bitmask = bitmsk(bits);        /* bitmask within a single unit */
  371.     qprec = bits2units(bits+1);
  372.     mp_setbit(remainder,(bits-RECIPMARGIN)-1);
  373.     /* rescale quotient to precision of divisor + RECIPMARGIN bits */
  374.     rescale(quotient,global_precision,qprec);
  375.     make_msbptr(quotient,qprec); 
  376.  
  377.     while (bits--)
  378.     {    mp_shift_left(remainder);
  379.         if (mp_compare(remainder,divisor) >= 0)
  380.         {    mp_sub(remainder,divisor);
  381.             stuff_bit(quotient,bitmask);
  382.         }
  383.         bump_bitsniffer(quotient,bitmask);
  384.     }
  385.     mp_init0(remainder);    /* burn sensitive data left on stack */
  386.     return(0);
  387. } /* mp_recip */
  388.  
  389.  
  390. int mp_div(register unitptr remainder,register unitptr quotient,
  391.     register unitptr dividend,register unitptr divisor)
  392.     /* Signed divide, either or both operands may be negative. */
  393. {    boolean dvdsign,dsign;
  394.     int status;
  395.     dvdsign = (mp_tstminus(dividend)!=0);
  396.     dsign = (mp_tstminus(divisor)!=0);
  397.     if (dvdsign) mp_neg(dividend);
  398.     if (dsign) mp_neg(divisor);
  399.     status = mp_udiv(remainder,quotient,dividend,divisor);
  400.     if (dvdsign) mp_neg(dividend);        /* repair caller's dividend */
  401.     if (dsign) mp_neg(divisor);        /* repair caller's divisor */
  402.     if (status<0) return(status);        /* divide error? */
  403.     if (dvdsign) mp_neg(remainder);
  404.     if (dvdsign ^ dsign) mp_neg(quotient);
  405.     return(status);
  406. } /* mp_div */
  407.  
  408.  
  409. word16 mp_shortdiv(register unitptr quotient,
  410.     register unitptr dividend,register word16 divisor)
  411. /*    This function does a fast divide and mod on a multiprecision dividend
  412.     using a short integer divisor returning a short integer remainder.
  413.     This is an unsigned divide.  It treats both operands as positive.
  414.     It is used mainly for faster printing of large numbers in base 10. 
  415. */
  416. {    int bits;
  417.     short dprec;
  418.     register unit bitmask;
  419.     register word16 remainder;
  420.     if (!divisor)    /* if divisor == 0 */
  421.         return(-1);    /* zero divisor means divide error */
  422.     remainder=0;
  423.     mp_init0(quotient);
  424.     /* normalize and compute number of bits in dividend first */
  425.     init_bitsniffer(dividend,bitmask,dprec,bits);
  426.     /* rescale quotient to same precision (dprec) as dividend */
  427.     rescale(quotient,global_precision,dprec);
  428.     make_msbptr(quotient,dprec); 
  429.  
  430.     while (bits--)
  431.     {    remainder <<= 1;
  432.         if (sniff_bit(dividend,bitmask))
  433.             remainder++;
  434.         if (remainder >= divisor)
  435.         {    remainder -= divisor;
  436.             stuff_bit(quotient,bitmask);
  437.         }
  438.         bump_2bitsniffers(dividend,quotient,bitmask); 
  439.     }    
  440.     return(remainder);
  441. } /* mp_shortdiv */
  442.  
  443.  
  444. int mp_mod(register unitptr remainder,
  445.     register unitptr dividend,register unitptr divisor)
  446.     /* Unsigned divide, treats both operands as positive. */
  447. {    int bits;
  448.     short dprec;
  449.     register unit bitmask;
  450.     if (testeq(divisor,0))
  451.         return(-1);    /* zero divisor means divide error */
  452.     mp_init0(remainder);
  453.     /* normalize and compute number of bits in dividend first */
  454.     init_bitsniffer(dividend,bitmask,dprec,bits);
  455.  
  456.     while (bits--)
  457.     {    mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
  458.         msub(remainder,divisor);
  459.         bump_bitsniffer(dividend,bitmask);
  460.     }    
  461.     return(0);
  462. } /* mp_mod */
  463.  
  464.  
  465. word16 mp_shortmod(register unitptr dividend,register word16 divisor)
  466. /*    This function does a fast mod operation on a multprecision dividend
  467.     using a short integer modulus returning a short integer remainder.
  468.     This is an unsigned divide.  It treats both operands as positive.
  469.     It is used mainly for fast sieve searches for large primes. 
  470. */
  471. {    int bits;
  472.     short dprec;
  473.     register unit bitmask;
  474.     register word16 remainder;
  475.     if (!divisor)    /* if divisor == 0 */
  476.         return(-1);    /* zero divisor means divide error */
  477.     remainder=0;
  478.     /* normalize and compute number of bits in dividend first */
  479.     init_bitsniffer(dividend,bitmask,dprec,bits);
  480.  
  481.     while (bits--)
  482.     {    remainder <<= 1;
  483.         if (sniff_bit(dividend,bitmask))
  484.             remainder++;
  485.         if (remainder >= divisor) remainder -= divisor;
  486.         bump_bitsniffer(dividend,bitmask);
  487.     }    
  488.     return(remainder);
  489. } /* mp_shortmod */
  490.  
  491.  
  492.  
  493. #ifdef COMB_MULT    /* use faster "comb" multiply algorithm */
  494.     /* We are skipping this code because it has a bug... */
  495.  
  496. int mp_mult(register unitptr prod,
  497.     register unitptr multiplicand, register unitptr multiplier)
  498.     /*    Computes multiprecision prod = multiplicand * multiplier */
  499. {    /*    Uses interleaved comb multiply algorithm.
  500.         This improved multiply more than twice as fast as a Russian 
  501.         peasant multiply, because it does a lot fewer shifts. 
  502.         Must have global_precision set to the size of the multiplicand 
  503.         plus UNITSIZE-1 SLOP_BITS.  Produces a product that is the sum 
  504.         of the lengths of the multiplier and multiplicand.
  505.  
  506.         BUG ALERT:  Unfortunately, this code has a bug.  It fails for 
  507.         some numbers.  One such example:
  508.         x=   59DE 60CE 2345 8091 A02B 2A1C DBC3 8BE5 
  509.         x*x= 59DE 60CE 2345 26B3 993B 67A5 2499 0B7D 
  510.              52C8 CDC7 AFB3 61C8 243C 741B
  511.         --which is obviously wrong.  The answer should be:
  512.         x*x= 1F8C 607B 5EA6 C061 2714 04A9 A0C6 A17A 
  513.              C9AB 6095 C62F 3756 3843 E4D0 3950 7AD9 
  514.         We'll have to fix this some day.  In the meantime, we'll 
  515.         just have the compiler skip over this code. 
  516.     */
  517.     int bits;
  518.     register unit bitmask;
  519.     unitptr product, mplier, temp;
  520.     short mprec,mprec2;
  521.     unit mplicand[MAX_UNIT_PRECISION];
  522.     mp_init(prod,0);        /* better clear full width--double precision */
  523.     if (testeq(multiplicand,0))
  524.         return(0);    /* zero multiplicand means zero product */
  525.  
  526.     mp_move(mplicand,multiplicand);    /* save it from damage */
  527.  
  528.     normalize(multiplier,mprec);
  529.     if (!mprec)
  530.         return(0);
  531.  
  532.     make_lsbptr(multiplier,mprec);
  533.     bitmask = 1;    /* start scan at LSB of multiplier */
  534.  
  535.     do    /* UNITSIZE times */
  536.     {    /* do for bits 0-15 */
  537.         product = prod;
  538.         mplier = multiplier;
  539.         mprec2 = mprec;
  540.         while (mprec2--)    /* do for each word in multiplier */
  541.         {
  542.             if (sniff_bit(mplier,bitmask))
  543.             {    if (mp_addc(product,multiplicand,0))    /* ripple carry */
  544.                 {    /* After 1st time thru, this is rarely encountered. */
  545.                     temp = msbptr(product,global_precision);
  546.                     pre_higherunit(temp);
  547.                     /* temp now points to LSU of carry region. */
  548.                     unmake_lsbptr(temp,global_precision);
  549.                     mp_inc(temp);
  550.                 }    /* ripple carry */
  551.             }
  552.             pre_higherunit(mplier);
  553.             pre_higherunit(product);
  554.         }
  555.         if (!(bitmask <<= 1))
  556.             break;
  557.         mp_shift_left(multiplicand);
  558.  
  559.     } while (TRUE);
  560.  
  561.     mp_move(multiplicand,mplicand);    /* recover */
  562.  
  563.     return(0);    /* normal return */    
  564. }    /* mp_mult */
  565.  
  566. #endif    /* COMB_MULT */
  567.  
  568.  
  569. /*    Because the "comb" multiply has a bug, we will use the slower
  570.     Russian peasant multiply instead.  Fortunately, the mp_mult 
  571.     function is not called from any time-critical code.
  572. */
  573.  
  574. int mp_mult(register unitptr prod,
  575.     register unitptr multiplicand,register unitptr multiplier)
  576.     /* Computes multiprecision prod = multiplicand * multiplier */
  577. {    /* Uses "Russian peasant" multiply algorithm. */
  578.     int bits;
  579.     register unit bitmask;
  580.     short mprec;
  581.     mp_init(prod,0);
  582.     if (testeq(multiplicand,0))
  583.         return(0);    /* zero multiplicand means zero product */
  584.     /* normalize and compute number of bits in multiplier first */
  585.     init_bitsniffer(multiplier,bitmask,mprec,bits);
  586.  
  587.     while (bits--)
  588.     {    mp_shift_left(prod);
  589.         if (sniff_bit(multiplier,bitmask))
  590.             mp_add(prod,multiplicand);
  591.         bump_bitsniffer(multiplier,bitmask);
  592.     }
  593.     return(0);    
  594. }    /* mp_mult */
  595.  
  596.  
  597.  
  598. /*    mp_modmult computes a multiprecision multiply combined with a 
  599.     modulo operation.  This is the most time-critical function in
  600.     this multiprecision arithmetic library for performing modulo
  601.     exponentiation.  We experimented with different versions of modmult,
  602.     depending on the machine architecture and performance requirements.
  603.     We will either use a Russian Peasant modmult (peasant_modmult), 
  604.     Charlie Merritt's modmult (merritt_modmult) or Jimmy Upton's
  605.     modmult (upton_modmult).  On machines with a hardware atomic 
  606.     multiply instruction, Upton's modmult is fastest, which depends
  607.     on an assembly subroutine to speed up the hardware multiply logic.
  608.     If the machine lacks a fast hardware multiply, Merritt's modmult
  609.     is preferred, which doesn't call any assembly multiply routine.
  610.     We use the alias names mp_modmult, stage_modulus, and modmult_burn 
  611.     for the corresponding true names, which depend on what flavor of 
  612.     modmult we are using.
  613.  
  614.     Before making the first call to mp_modmult, you must set up the 
  615.     modulus-dependant precomputated tables by calling stage_modulus.
  616.     After making all the calls to mp_modmult, you call modmult_burn to 
  617.     erase the tables created by stage_modulus that were left in memory.
  618. */
  619.  
  620. #ifdef COUNTMULTS
  621. /* "number of modmults" counters, used for performance studies. */
  622. static unsigned int tally_modmults = 0;
  623. static unsigned int tally_modsquares = 0;
  624. #endif    /* COUNTMULTS */
  625.  
  626.  
  627. #ifdef PEASANT
  628. /* Conventional Russian peasant multiply with modulo algorithm. */
  629.  
  630. static unitptr modulus = 0;    /* used only by mp_modmult */
  631.  
  632. int stage_peasant_modulus(unitptr n)
  633. /*    Must pass modulus to stage_modulus before calling modmult.
  634.     Assumes that global_precision has already been adjusted to the
  635.     size of the modulus, plus SLOP_BITS.
  636. */
  637. {    /* For this simple version of modmult, just copy unit pointer. */
  638.     modulus = n;
  639.     return(0);    /* normal return */
  640. }    /* stage_peasant_modulus */
  641.  
  642.  
  643. int peasant_modmult(register unitptr prod,
  644.     unitptr multiplicand,register unitptr multiplier)
  645. {    /*    "Russian peasant" multiply algorithm, combined with a modulo 
  646.         operation.  This is a simple naive replacement for Merritt's 
  647.         faster modmult algorithm.  References global unitptr "modulus".
  648.         Computes:  prod = (multiplicand*multiplier) mod modulus
  649.         WARNING: All the arguments must be less than the modulus!
  650.     */
  651.     int bits;
  652.     register unit bitmask;
  653.     short mprec;
  654.     mp_init(prod,0);
  655. /*    if (testeq(multiplicand,0))
  656.         return(0); */    /* zero multiplicand means zero product */
  657.     /* normalize and compute number of bits in multiplier first */
  658.     init_bitsniffer(multiplier,bitmask,mprec,bits);
  659.  
  660.     while (bits--)
  661.     {    mp_shift_left(prod);
  662.         msub(prod,modulus);    /* turns mult into modmult */
  663.         if (sniff_bit(multiplier,bitmask))
  664.         {    mp_add(prod,multiplicand);
  665.             msub(prod,modulus);    /* turns mult into modmult */
  666.         }
  667.         bump_bitsniffer(multiplier,bitmask);
  668.     }
  669.     return(0);
  670. }    /* peasant_modmult */
  671.  
  672.  
  673. /*    If we are using a version of mp_modmult that uses internal tables 
  674.     in memory, we have to call modmult_burn() at the end of mp_modexp.
  675.     This is so that no sensitive data is left in memory after the program 
  676.     exits.  The Russian peasant method doesn't use any such tables.
  677. */
  678. static void peasant_burn(void)
  679. /*    Alias for modmult_burn, called only from mp_modexp().  Destroys
  680.     internal modmult tables.  This version does nothing because no 
  681.     tables are used by the Russian peasant modmult. */
  682. { }    /* peasant_burn */
  683.  
  684. #endif    /* PEASANT */
  685.  
  686.  
  687. #ifdef MERRITT
  688. /*=========================================================================*/
  689. /*
  690.     This is Charlie Merritt's MODMULT algorithm, implemented in C by PRZ.
  691.     Also refined by Zhahai Stewart to reduce the number of subtracts.
  692.     It performs a multiply combined with a modulo operation, without 
  693.     going into "double precision".  It is faster than the Russian peasant
  694.     method,    and still works well on machines that lack a fast hardware 
  695.     multiply instruction.  
  696. */
  697.  
  698. /*    The following support functions, macros, and data structures
  699.     are used only by Merritt's modmult algorithm... */
  700.  
  701. static void mp_lshift_unit(register unitptr r1)
  702. /*    Shift r1 1 whole unit to the left.  Used only by modmult function. */
  703. {    register short precision;
  704.     register unitptr r2;
  705.     precision = global_precision;
  706.     make_msbptr(r1,precision);
  707.     r2 = r1;
  708.     while (--precision)
  709.         *post_lowerunit(r1) = *pre_lowerunit(r2);
  710.     *r1 = 0;
  711. }    /* mp_lshift_unit */
  712.  
  713.  
  714. /* moduli_buf contains shifted images of the modulus, set by stage_modulus */
  715. static unit moduli_buf[UNITSIZE][MAX_UNIT_PRECISION] = {0};
  716. static unitptr moduli[UNITSIZE+1] = /* contains pointers into moduli_buf */
  717. {    0
  718.     ,&moduli_buf[ 0][0], &moduli_buf[ 1][0], &moduli_buf[ 2][0], &moduli_buf[ 3][0], 
  719.      &moduli_buf[ 4][0], &moduli_buf[ 5][0], &moduli_buf[ 6][0], &moduli_buf[ 7][0]
  720. #ifndef UNIT8
  721.     ,&moduli_buf[ 8][0], &moduli_buf[ 9][0], &moduli_buf[10][0], &moduli_buf[11][0], 
  722.      &moduli_buf[12][0], &moduli_buf[13][0], &moduli_buf[14][0], &moduli_buf[15][0] 
  723. #ifndef UNIT16    /* and not UNIT8 */
  724.     ,&moduli_buf[16][0], &moduli_buf[17][0], &moduli_buf[18][0], &moduli_buf[19][0], 
  725.      &moduli_buf[20][0], &moduli_buf[21][0], &moduli_buf[22][0], &moduli_buf[23][0], 
  726.      &moduli_buf[24][0], &moduli_buf[25][0], &moduli_buf[26][0], &moduli_buf[27][0], 
  727.      &moduli_buf[28][0], &moduli_buf[29][0], &moduli_buf[30][0], &moduli_buf[31][0]
  728. #endif    /* UNIT16 and UNIT8 not defined */
  729. #endif    /* UNIT8 not defined */
  730. };
  731.  
  732. /*    To optimize msubs, need following 2 unit arrays, each filled
  733.     with the most significant unit and next-to-most significant unit
  734.     of the preshifted images of the modulus. */
  735. static unit msu_moduli[UNITSIZE+1] = {0}; /* most signif. unit */
  736. static unit nmsu_moduli[UNITSIZE+1] = {0}; /* next-most signif. unit */
  737.  
  738. /*    mpdbuf contains preshifted images of the multiplicand, mod n.
  739.      It is used only by mp_modmult.  It could be staticly declared
  740.     inside of mp_modmult, but we put it outside mp_modmult so that 
  741.     it can be wiped clean by modmult_burn(), which is called at the
  742.     end of mp_modexp.  This is so that no sensitive data is left in 
  743.     memory after the program exits.
  744. */
  745. static unit mpdbuf[UNITSIZE-1][MAX_UNIT_PRECISION] = {0};
  746.  
  747.  
  748. static void stage_mp_images(unitptr images[UNITSIZE],unitptr r)
  749. /*    Computes UNITSIZE images of r, each shifted left 1 more bit.
  750.     Used only by modmult function. 
  751. */
  752. {    short int i;
  753.     images[0] = r;    /* no need to move the first image, just copy ptr */
  754.     for (i=1; i<UNITSIZE; i++)
  755.     {    mp_move(images[i],images[i-1]);
  756.         mp_shift_left(images[i]);
  757.     }
  758. }    /* stage_mp_images */
  759.  
  760.  
  761. int stage_merritt_modulus(unitptr n)
  762. /*    Computes UNITSIZE+1 images of modulus, each shifted left 1 more bit.
  763.     Before calling mp_modmult, you must first stage the moduli images by
  764.     calling stage_modulus.  n is the pointer to the modulus.
  765.     Assumes that global_precision has already been adjusted to the
  766.     size of the modulus, plus SLOP_BITS.
  767. */
  768. {    short int i;
  769.     unitptr msu;    /* ptr to most significant unit, for faster msubs */
  770.     moduli[0] = n;    /* no need to move the first image, just copy ptr */
  771.  
  772.     /* used by optimized msubs macro... */
  773.     msu = msbptr(n,global_precision);    /* needed by msubs */
  774.     msu_moduli[0] = *post_lowerunit(msu);
  775.     nmsu_moduli[0] = *msu;
  776.  
  777.     for (i=1; i<UNITSIZE+1; i++)
  778.     {    mp_move(moduli[i],moduli[i-1]);
  779.         mp_shift_left(moduli[i]);
  780.  
  781.         /* used by optimized msubs macro... */
  782.         msu = msbptr(moduli[i],global_precision);    /* needed by msubs */
  783.         msu_moduli[i] = *post_lowerunit(msu);    /* for faster msubs */
  784.         nmsu_moduli[i] = *msu;
  785.     }
  786.     return(0);    /* normal return */
  787. }    /* stage_merritt_modulus */
  788.  
  789.  
  790. /* The following macros, sniffadd and msubs, are used by modmult... */
  791. #define sniffadd(i) if (*multiplier & power_of_2(i))  mp_add(prod,mpd[i])
  792.  
  793. /* Unoptimized msubs macro (msubs0) follows... */
  794. /* #define msubs0(i) msub(prod,moduli[i]) 
  795. */
  796.  
  797. /*    To optimize msubs, compare the most significant units of the 
  798.     product and the shifted modulus before deciding to call the full 
  799.     mp_compare routine.  Better still, compare the upper two units
  800.     before deciding to call mp_compare.
  801.     Optimization of msubs relies heavily on standard C left-to-right 
  802.     parsimonious evaluation of logical expressions. 
  803. */
  804.  
  805. /* Partially-optimized msubs macro (msubs1) follows... */
  806. /* #define msubs1(i) if ( \
  807.   ((p_m = (*msu_prod-msu_moduli[i])) >= 0) && \
  808.   (p_m || (mp_compare(prod,moduli[i]) >= 0) ) \
  809.   ) mp_sub(prod,moduli[i])
  810. */
  811.  
  812. /* Fully-optimized msubs macro (msubs2) follows... */
  813. #define msubs(i) if (((p_m = *msu_prod-msu_moduli[i])>0) || ( \
  814.  (p_m==0) && ( (*nmsu_prod>nmsu_moduli[i]) || ( \
  815.  (*nmsu_prod==nmsu_moduli[i]) && ((mp_compare(prod,moduli[i]) >= 0)) ))) ) \
  816.  mp_sub(prod,moduli[i])
  817.  
  818.  
  819. int merritt_modmult(register unitptr prod,
  820.     unitptr multiplicand,register unitptr multiplier)
  821.     /*    Performs combined multiply/modulo operation.  
  822.         Computes:  prod = (multiplicand*multiplier) mod modulus
  823.         WARNING: All the arguments must be less than the modulus!
  824.         Assumes the modulus has been predefined by first calling 
  825.         stage_modulus.
  826.     */
  827. {
  828.     /* p_m, msu_prod, and nmsu_prod are used by the optimized msubs macro...*/
  829.     register signedunit p_m;
  830.     register unitptr msu_prod;    /* ptr to most significant unit of product */
  831.     register unitptr nmsu_prod;    /* next-most signif. unit of product */
  832.     short mprec;        /* precision of multiplier, in units */
  833.     /*    Array mpd contains a list of pointers to preshifted images of 
  834.         the multiplicand: */
  835.     static unitptr mpd[UNITSIZE] =
  836.     {     0,              &mpdbuf[ 0][0], &mpdbuf[ 1][0], &mpdbuf[ 2][0],
  837.          &mpdbuf[ 3][0], &mpdbuf[ 4][0], &mpdbuf[ 5][0], &mpdbuf[ 6][0]
  838. #ifndef UNIT8
  839.         ,&mpdbuf[ 7][0], &mpdbuf[ 8][0], &mpdbuf[ 9][0], &mpdbuf[10][0],
  840.          &mpdbuf[11][0], &mpdbuf[12][0], &mpdbuf[13][0], &mpdbuf[14][0]
  841. #ifndef UNIT16    /* and not UNIT8 */
  842.         ,&mpdbuf[15][0], &mpdbuf[16][0], &mpdbuf[17][0], &mpdbuf[18][0],
  843.          &mpdbuf[19][0], &mpdbuf[20][0], &mpdbuf[21][0], &mpdbuf[22][0],
  844.          &mpdbuf[23][0], &mpdbuf[24][0], &mpdbuf[25][0], &mpdbuf[26][0],
  845.          &mpdbuf[27][0], &mpdbuf[28][0], &mpdbuf[29][0], &mpdbuf[30][0]
  846. #endif    /* UNIT16 and UNIT8 not defined */
  847. #endif    /* UNIT8 not defined */
  848.     };
  849.  
  850.     /* Compute preshifted images of multiplicand, mod n: */
  851.     stage_mp_images(mpd,multiplicand);
  852.  
  853.     /* To optimize msubs, set up msu_prod and nmsu_prod: */
  854.     msu_prod = msbptr(prod,global_precision); /* Get ptr to MSU of prod */
  855.     nmsu_prod = msu_prod;
  856.     post_lowerunit(nmsu_prod); /* Get next-MSU of prod */
  857.  
  858.     /*    To understand this algorithm, it would be helpful to first 
  859.         study the conventional Russian peasant modmult algorithm.
  860.         This one does about the same thing as Russian peasant, but
  861.         is organized differently to save some steps.  It loops 
  862.         through the multiplier a word (unit) at a time, instead of 
  863.         a bit at a time.  It word-shifts the product instead of 
  864.         bit-shifting it, so it should be faster.  It also does about 
  865.         half as many subtracts as Russian peasant.
  866.     */
  867.  
  868.     mp_init(prod,0);    /* Initialize product to 0. */
  869.  
  870.     /*    The way mp_modmult is actually used in cryptographic 
  871.         applications, there will NEVER be a zero multiplier or 
  872.         multiplicand.  So there is no need to optimize for that 
  873.         condition.
  874.     */
  875. /*    if (testeq(multiplicand,0))
  876.         return(0); */    /* zero multiplicand means zero product */
  877.     /* Normalize and compute number of units in multiplier first: */
  878.     normalize(multiplier,mprec);
  879.     if (mprec==0)    /* if precision of multiplier is 0 */
  880.         return(0);    /* zero multiplier means zero product */
  881.     make_msbptr(multiplier,mprec); /* start at MSU of multiplier */
  882.  
  883.     while (mprec--)    /* Loop for the number of units in the multiplier */
  884.     {
  885.         /*    Shift the product one whole unit to the left.
  886.             This is part of the multiply phase of modmult. 
  887.         */
  888.  
  889.         mp_lshift_unit(prod);
  890.  
  891.         /*    Sniff each bit in the current unit of the multiplier, 
  892.             and conditionally add the the corresponding preshifted 
  893.             image of the multiplicand to the product.
  894.             This is also part of the multiply phase of modmult.
  895.  
  896.             The following loop is unrolled for speed, using macros...
  897.  
  898.         for (i=UNITSIZE-1; i>=0; i--)
  899.            if (*multiplier & power_of_2(i)) 
  900.                 mp_add(prod,mpd[i]);
  901.         */
  902. #ifndef UNIT8
  903. #ifndef UNIT16    /* and not UNIT8 */
  904.         sniffadd(31); 
  905.         sniffadd(30); 
  906.         sniffadd(29); 
  907.         sniffadd(28);
  908.         sniffadd(27); 
  909.         sniffadd(26); 
  910.         sniffadd(25); 
  911.         sniffadd(24);
  912.         sniffadd(23); 
  913.         sniffadd(22); 
  914.         sniffadd(21); 
  915.         sniffadd(20);
  916.         sniffadd(19); 
  917.         sniffadd(18); 
  918.         sniffadd(17); 
  919.         sniffadd(16);
  920. #endif    /* not UNIT16 and not UNIT8 */
  921.         sniffadd(15); 
  922.         sniffadd(14); 
  923.         sniffadd(13); 
  924.         sniffadd(12);
  925.         sniffadd(11); 
  926.         sniffadd(10); 
  927.         sniffadd(9); 
  928.         sniffadd(8);
  929. #endif    /* not UNIT8 */
  930.         sniffadd(7); 
  931.         sniffadd(6); 
  932.         sniffadd(5); 
  933.         sniffadd(4);
  934.         sniffadd(3); 
  935.         sniffadd(2); 
  936.         sniffadd(1); 
  937.         sniffadd(0);
  938.  
  939.         /*    The product may have grown by as many as UNITSIZE+1 
  940.             bits.  That's why we have global_precision set to the
  941.             size of the modulus plus UNITSIZE+1 slop bits.  
  942.             Now reduce the product back down by conditionally 
  943.             subtracting    the UNITSIZE+1 preshifted images of the 
  944.             modulus.  This is the modulo reduction phase of modmult.
  945.  
  946.             The following loop is unrolled for speed, using macros...
  947.  
  948.         for (i=UNITSIZE; i>=0; i--)
  949.             if (mp_compare(prod,moduli[i]) >= 0) 
  950.                 mp_sub(prod,moduli[i]); 
  951.         */
  952. #ifndef UNIT8
  953. #ifndef UNIT16    /* and not UNIT8 */
  954.         msubs(32); 
  955.         msubs(31); 
  956.         msubs(30); 
  957.         msubs(29); 
  958.         msubs(28);
  959.         msubs(27); 
  960.         msubs(26); 
  961.         msubs(25); 
  962.         msubs(24);
  963.         msubs(23); 
  964.         msubs(22); 
  965.         msubs(21); 
  966.         msubs(20);
  967.         msubs(19); 
  968.         msubs(18); 
  969.         msubs(17); 
  970. #endif    /* not UNIT16 and not UNIT8 */
  971.         msubs(16);
  972.         msubs(15); 
  973.         msubs(14); 
  974.         msubs(13); 
  975.         msubs(12);
  976.         msubs(11); 
  977.         msubs(10); 
  978.         msubs(9); 
  979. #endif    /* not UNIT8 */
  980.         msubs(8);
  981.         msubs(7); 
  982.         msubs(6); 
  983.         msubs(5); 
  984.         msubs(4);
  985.         msubs(3); 
  986.         msubs(2); 
  987.         msubs(1); 
  988.         msubs(0);
  989.  
  990.         /* Bump pointer to next lower unit of multiplier: */
  991.         post_lowerunit(multiplier);
  992.  
  993.     }    /* Loop for the number of units in the multiplier */
  994.  
  995.     return(0);    /* normal return */
  996.  
  997. }    /* merritt_modmult */
  998.  
  999.  
  1000. #undef msubs
  1001. #undef sniffadd
  1002.  
  1003.  
  1004. /*    Merritt's mp_modmult function leaves some internal tables in memory,
  1005.     so we have to call modmult_burn() at the end of mp_modexp.  
  1006.     This is so that no cryptographically sensitive data is left in memory 
  1007.     after the program exits.
  1008. */
  1009. static void merritt_burn(void)
  1010. /*    Alias for modmult_burn, merritt_burn() is called only by mp_modexp. */
  1011. {    unitfill0(&(mpdbuf[0][0]),(UNITSIZE-1)*MAX_UNIT_PRECISION);
  1012.     unitfill0(&(moduli_buf[0][0]),(UNITSIZE)*MAX_UNIT_PRECISION);
  1013.     unitfill0(msu_moduli,UNITSIZE+1);
  1014.     unitfill0(nmsu_moduli,UNITSIZE+1);
  1015. } /* merritt_burn() */
  1016.  
  1017. /******* end of Merritt's MODMULT stuff. *******/
  1018. /*=========================================================================*/
  1019. #endif    /* MERRITT */
  1020.  
  1021.  
  1022. #ifdef UPTON    /* using Jimmy Upton's modmult algorithm */
  1023.  
  1024. /*    Jimmy Upton's multiprecision modmult algorithm in C.
  1025.     Performs a multiply combined with a modulo operation.  
  1026.  
  1027.     The following support functions and data structures
  1028.     are used only by Upton's modmult algorithm...
  1029. */
  1030.  
  1031. /*    The MULTUNIT word is the biggest word size for the native atomic 
  1032.     multiply instruction.  It may or may not be smaller than UNITSIZE. 
  1033.     Many CPU's have 16-by-16-bit multiplies, yielding a 32-bit product.
  1034.     In that case, MULTUINIT should be a 16-bit word, even if the rest of
  1035.     the multiprecision arithmetic functions assume a 32-bit UNIT word.
  1036. */
  1037. typedef unsigned short MULTUNIT;
  1038. #define    MULTUNITSIZE    (8*sizeof(MULTUNIT))    /* size in bits */
  1039.  
  1040. static short unit_prec;    /* global_precision expressed in MULTUNITs */
  1041.  
  1042.  
  1043. /*    Define MPORTABLE if there is no assembly version of the mp_smul
  1044.     function available.  An assembly version is much faster.
  1045.     Note that since the SPARC CPU has no hardware integer multiply 
  1046.     instruction, there is not that much advantage in having an 
  1047.     assembly version of mp_smul on that machine.  It might be faster
  1048.     to use Merritt's modmult instead of Upton's modmult on the SPARC.
  1049. */
  1050. #ifdef MSDOS    /* we do have an assembly version available on the 8086 */
  1051. #undef MPORTABLE
  1052. #endif
  1053.  
  1054. #ifdef MPORTABLE  /* use slow portable C version instead of assembly */
  1055.  
  1056. /* 
  1057.     Multiply the single-word multiplier times the multiprecision integer 
  1058.     in multiplicand, accumulating result in prod.  The resulting 
  1059.     multiprecision prod will be 1 word longer than the multiplicand.   
  1060.     multiplicand is unit_prec words long.  We add into prod, so caller 
  1061.     should zero it out first.  For best results, this time-critical 
  1062.     function should be implemented in assembly.
  1063.     NOTE:  Unlike other functions in the multiprecision arithmetic 
  1064.     library, both multiplicand and prod are pointing at the LSB, 
  1065.     regardless of byte order of the machine.  On an 80x86, this makes 
  1066.     no difference.  But if this assembly function is implemented
  1067.     on a 680x0, it becomes important.
  1068. */
  1069. void mp_smul (MULTUNIT *prod, MULTUNIT *multiplicand, MULTUNIT multiplier)
  1070. {
  1071.     short i;
  1072.     unsigned long p, carry;
  1073.  
  1074.     carry = 0;
  1075.     for (i=0; i<unit_prec; ++i)
  1076.     {    p = (unsigned long)multiplier * *post_higherunit(multiplicand);
  1077.         p += *prod + carry;
  1078.         *post_higherunit(prod) = (MULTUNIT)p;
  1079.         carry = p >> MULTUNITSIZE;
  1080.     }
  1081.     /* We know that the high-order word of prod will always be 0 */
  1082.     *prod = (MULTUNIT)carry;    /* copy carry to empty high word of prod */
  1083. }    /* mp_smul */
  1084.  
  1085. #else    /* not MPORTABLE-- use assembly routine */
  1086.  
  1087. #define mp_smul P_SMUL    /* use external assembly routine */
  1088.  
  1089. #endif /* MPORTABLE */
  1090.  
  1091. #ifdef VMS
  1092.  
  1093. #define mp_dmul p_dmul /* use external assembly routine */
  1094.  
  1095. void mp_dmul (unitptr prod, unitptr multiplicand, unitptr multiplier);
  1096.  
  1097. #else /* VMS */
  1098.  
  1099. /*    mp_dmul is a double-precision multiply multiplicand times multiplier, 
  1100.     result into prod.  prod must be pointing at a "double precision" 
  1101.     buffer.  E.g. If global_precision is 10 words, prod must be 
  1102.     pointing at a 20-word buffer.
  1103. */
  1104. static void mp_dmul (unitptr prod, unitptr multiplicand, unitptr multiplier)
  1105. {
  1106.     register    int i;
  1107.     register    MULTUNIT *p_multiplicand, *p_multiplier;
  1108.     register    MULTUNIT *prodp;
  1109.  
  1110.     unitfill0(prod,global_precision*2);    /* Pre-zero prod */
  1111.     /* Calculate precision in units of MULTUNIT */
  1112.     unit_prec = global_precision * UNITSIZE / MULTUNITSIZE;
  1113.     p_multiplicand = (MULTUNIT *)multiplicand;
  1114.     p_multiplier = (MULTUNIT *)multiplier;
  1115.     prodp = (MULTUNIT *)prod;
  1116.     make_lsbptr(p_multiplicand,unit_prec);
  1117.     make_lsbptr(p_multiplier,unit_prec);
  1118.     make_lsbptr(prodp,unit_prec*2);
  1119.     /* Multiply multiplicand by each word in multiplier, accumulating prod: */
  1120.     for (i=0; i<unit_prec; ++i)
  1121.         mp_smul (post_higherunit(prodp), 
  1122.             p_multiplicand, *post_higherunit(p_multiplier));
  1123. }    /* mp_dmul */
  1124.  
  1125. #endif /* VMS */
  1126.  
  1127. /*    These scratchpad arrays are used only by upton_modmult (mp_modmult).
  1128.      Some of them could be staticly declared inside of mp_modmult, but we 
  1129.     put them outside mp_modmult so that they can be wiped clean by 
  1130.     modmult_burn(), which is called at the end of mp_modexp.  This is 
  1131.     so that no sensitive data is left in memory after the program exits.
  1132. */
  1133.  
  1134. #ifdef VAXC
  1135. /*
  1136.  * Operations on arrays are sped up on a VAX if the data is quadword (64-bit)
  1137.  * aligned. This is because many VAXen use a minimum of 64-bit data paths to
  1138.  * cache and memory. This technique is VAX C private, hence the conditional.
  1139.  */
  1140. static unit _align( quadword ) modulus[MAX_UNIT_PRECISION];
  1141. static unit _align( quadword ) reciprocal[MAX_UNIT_PRECISION];
  1142. static unit _align( quadword ) dhi[MAX_UNIT_PRECISION];
  1143. static unit _align( quadword ) d_data[MAX_UNIT_PRECISION*2];     /* double-wide register */
  1144. static unit _align( quadword ) e_data[MAX_UNIT_PRECISION*2];     /* double-wide register */
  1145. static unit _align( quadword ) f_data[MAX_UNIT_PRECISION*2];     /* double-wide register */
  1146.  
  1147. #else /* VAXC */
  1148.  
  1149. static unit modulus[MAX_UNIT_PRECISION];
  1150. static unit reciprocal[MAX_UNIT_PRECISION];
  1151. static unit dhi[MAX_UNIT_PRECISION];
  1152. static unit d_data[MAX_UNIT_PRECISION*2];    /* double-wide register */
  1153. static unit e_data[MAX_UNIT_PRECISION*2];    /* double-wide register */
  1154. static unit f_data[MAX_UNIT_PRECISION*2];    /* double-wide register */
  1155.  
  1156. #endif /* VAXC */
  1157.  
  1158. static short nbits;
  1159. static short nbitsDivUNITSIZE;
  1160. static short nbitsModUNITSIZE;
  1161.  
  1162. /*    stage_upton_modulus() is aliased to stage_modulus().
  1163.     Prepare for an Upton modmult.  Calculate the reciprocal of modulus,
  1164.     and save both.  Note that reciprocal will have one more bit than 
  1165.     modulus.
  1166.     Assumes that global_precision has already been adjusted to the
  1167.     size of the modulus, plus SLOP_BITS.
  1168. */
  1169. int stage_upton_modulus(unitptr n)
  1170. {
  1171.     mp_move(modulus, n);
  1172.     mp_recip(reciprocal, modulus);
  1173.     nbits = countbits(modulus);
  1174.     nbitsDivUNITSIZE = nbits / UNITSIZE;
  1175.     nbitsModUNITSIZE = nbits % UNITSIZE;
  1176.     return(0);    /* normal return */
  1177. }    /* stage_upton_modulus */
  1178.  
  1179.  
  1180. /*     Upton's algorithm performs a multiply combined with a modulo operation.  
  1181.     Computes:  prod = (multiplicand*multiplier) mod modulus
  1182.     WARNING: All the arguments must be less than the modulus!
  1183.     References global unitptr modulus and reciprocal.
  1184.     The reciprocal of modulus is 1 bit longer than the modulus.  
  1185.     upton_modmult() is aliased to mp_modmult().
  1186. */
  1187. int upton_modmult (unitptr prod, unitptr multiplicand, unitptr multiplier)
  1188. {
  1189.     unitptr d = d_data;
  1190.     unitptr d1 = d_data;
  1191.     unitptr e = e_data;
  1192.     unitptr f = f_data;
  1193.     short orig_precision, dprec;
  1194.     short i;
  1195.  
  1196.     orig_precision = global_precision;
  1197.     mp_dmul (d,multiplicand,multiplier);
  1198.     
  1199.     /* Throw off low nbits of d */
  1200. #ifndef HIGHFIRST
  1201.     d1 = d + nbitsDivUNITSIZE;
  1202. #else
  1203.     d1 = d + orig_precision - nbitsDivUNITSIZE;
  1204. #endif
  1205.     mp_move(dhi, d1);    /* Don't screw up d, we need it later */
  1206.     mp_shift_right_bits(dhi, nbitsModUNITSIZE);
  1207.  
  1208.     mp_dmul (e,dhi,reciprocal);    /* Note - reciprocal has nbits+1 bits */
  1209.  
  1210. #ifndef HIGHFIRST
  1211.     e += nbitsDivUNITSIZE;
  1212. #else
  1213.     e += orig_precision - nbitsDivUNITSIZE;
  1214. #endif
  1215.     mp_shift_right_bits(e, nbitsModUNITSIZE);
  1216.  
  1217.     mp_dmul (f,e,modulus);
  1218.  
  1219.     /* Now for the only double-precision call to mpilib: */
  1220.     set_precision (orig_precision * 2);
  1221.     mp_sub (d,f);
  1222.  
  1223.     /* d's precision should be <= orig_precision */
  1224.     rescale (d, orig_precision*2, orig_precision);
  1225.     set_precision (orig_precision);
  1226.  
  1227.     /* Should never have to do this final subtract more than twice: */
  1228.     while (mp_compare(d,modulus) > 0)
  1229.         mp_sub (d,modulus);
  1230.  
  1231.     mp_move(prod,d);
  1232.     return(0);    /* normal return */
  1233. }    /* upton_modmult */
  1234.  
  1235.  
  1236. /*    Upton's mp_modmult function leaves some internal arrays in memory,
  1237.     so we have to call modmult_burn() at the end of mp_modexp.  
  1238.     This is so that no cryptographically sensitive data is left in memory 
  1239.     after the program exits.
  1240.     upton_burn() is aliased to modmult_burn().
  1241. */
  1242. void upton_burn(void)
  1243. {
  1244.     unitfill0(modulus,MAX_UNIT_PRECISION);
  1245.     unitfill0(reciprocal,MAX_UNIT_PRECISION);
  1246.     unitfill0(dhi,MAX_UNIT_PRECISION);
  1247.     unitfill0(d_data,MAX_UNIT_PRECISION*2);
  1248.     unitfill0(e_data,MAX_UNIT_PRECISION*2);
  1249.     unitfill0(f_data,MAX_UNIT_PRECISION*2);
  1250.     nbits = nbitsDivUNITSIZE = nbitsModUNITSIZE = 0;
  1251. }    /* upton_burn */
  1252.  
  1253. #endif    /* UPTON */
  1254.  
  1255.  
  1256. int countbits(unitptr r)
  1257.     /* Returns number of significant bits in r */
  1258. {    int bits;
  1259.     short prec;
  1260.     register unit bitmask;
  1261.     init_bitsniffer(r,bitmask,prec,bits);
  1262.     return(bits);
  1263. } /* countbits */
  1264.  
  1265.  
  1266. char *copyright_notice(void)
  1267. /* force linker to include copyright notice in the executable object image. */
  1268. { return ("(c)1986 Philip Zimmermann"); } /* copyright_notice */
  1269.  
  1270.  
  1271. int mp_modexp(register unitptr expout,register unitptr expin,
  1272.     register unitptr exponent,register unitptr modulus)
  1273. {    /*    Russian peasant combined exponentiation/modulo algorithm.
  1274.         Calls modmult instead of mult. 
  1275.         Computes:  expout = (expin**exponent) mod modulus
  1276.         WARNING: All the arguments must be less than the modulus!
  1277.     */
  1278.     int bits;
  1279.     short oldprecision;
  1280.     register unit bitmask;
  1281.     unit product[MAX_UNIT_PRECISION];
  1282.     short eprec;
  1283.  
  1284. #ifdef COUNTMULTS
  1285.     tally_modmults = 0;    /* clear "number of modmults" counter */
  1286.     tally_modsquares = 0;    /* clear "number of modsquares" counter */
  1287. #endif    /* COUNTMULTS */
  1288.     mp_init(expout,1);
  1289.     if (testeq(exponent,0))
  1290.     {    if (testeq(expin,0))
  1291.             return(-1); /* 0 to the 0th power means return error */
  1292.         return(0);    /* otherwise, zero exponent means expout is 1 */
  1293.     }
  1294.     if (testeq(modulus,0))
  1295.         return(-2);        /* zero modulus means error */
  1296. #if SLOP_BITS > 0    /* if there's room for sign bits */
  1297.     if (mp_tstminus(modulus))
  1298.         return(-2);        /* negative modulus means error */
  1299. #endif    /* SLOP_BITS > 0 */
  1300.     if (mp_compare(expin,modulus) >= 0)
  1301.         return(-3); /* if expin >= modulus, return error */
  1302.     if (mp_compare(exponent,modulus) >= 0)
  1303.         return(-4); /* if exponent >= modulus, return error */
  1304.  
  1305.     oldprecision = global_precision;    /* save global_precision */
  1306.     /* set smallest optimum precision for this modulus */
  1307.     set_precision(bits2units(countbits(modulus)+SLOP_BITS));
  1308.     /* rescale all these registers to global_precision we just defined */
  1309.     rescale(modulus,oldprecision,global_precision);
  1310.     rescale(expin,oldprecision,global_precision);
  1311.     rescale(exponent,oldprecision,global_precision);
  1312.     rescale(expout,oldprecision,global_precision);
  1313.  
  1314.     if (stage_modulus(modulus))
  1315.     {    set_precision(oldprecision);    /* restore original precision */
  1316.         return(-5);        /* unstageable modulus (STEWART algorithm) */
  1317.     }
  1318.  
  1319.     /* normalize and compute number of bits in exponent first */
  1320.     init_bitsniffer(exponent,bitmask,eprec,bits);
  1321.  
  1322.     /* We can "optimize out" the first modsquare and modmult: */
  1323.     bits--;        /* We know for sure at this point that bits>0 */
  1324.     mp_move(expout,expin);        /*  expout = (1*1)*expin; */
  1325.     bump_bitsniffer(exponent,bitmask);
  1326.     
  1327.     while (bits--)
  1328.     {
  1329.         poll_for_break(); /* polls keyboard, allows ctrl-C to abort program */
  1330. #ifdef COUNTMULTS
  1331.         tally_modsquares++;    /* bump "number of modsquares" counter */
  1332. #endif    /* COUNTMULTS */
  1333.         mp_modsquare(product,expout);
  1334.         mp_move(expout,product);
  1335.         if (sniff_bit(exponent,bitmask))
  1336.         {    mp_modmult(product,expout,expin);
  1337.             mp_move(expout,product);
  1338. #ifdef COUNTMULTS
  1339.             tally_modmults++;    /* bump "number of modmults" counter */
  1340. #endif    /* COUNTMULTS */
  1341.         }
  1342.         bump_bitsniffer(exponent,bitmask);
  1343.     }    /* while bits-- */
  1344.     mp_burn(product);    /* burn the evidence on the stack */
  1345.     modmult_burn(); /* ask mp_modmult to also burn its own evidence */
  1346.  
  1347. #ifdef COUNTMULTS    /* diagnostic analysis */
  1348.     {    long atomic_mults;
  1349.         unsigned int unitcount,totalmults;
  1350.         unitcount = bits2units(countbits(modulus));
  1351.         /* calculation assumes modsquare takes as long as a modmult: */
  1352.         atomic_mults = (long) tally_modmults * (unitcount * unitcount);
  1353.         atomic_mults += (long) tally_modsquares * (unitcount * unitcount);
  1354.         printf("%ld atomic mults for ",atomic_mults);
  1355.         printf("%d+%d = %d modsqr+modmlt, at %d bits, %d words.\n",
  1356.             tally_modsquares,tally_modmults,
  1357.             tally_modsquares+tally_modmults,
  1358.             countbits(modulus), unitcount);
  1359.     }
  1360. #endif    /* COUNTMULTS */
  1361.  
  1362.     set_precision(oldprecision);    /* restore original precision */
  1363.  
  1364.     /* Do an explicit reference to the copyright notice so that the linker
  1365.        will be forced to include it in the executable object image... */
  1366.     copyright_notice();    /* has no real effect at run time */
  1367.  
  1368.     return(0);        /* normal return */
  1369. }    /* mp_modexp */
  1370.  
  1371.  
  1372.  
  1373.  
  1374. /*********************************************************************
  1375.     RSA-specific routines follow.  These are the only functions that 
  1376.     are specific to the RSA public key cryptosystem.  The other
  1377.     multiprecision integer math functions may be used for non-RSA
  1378.     applications.  
  1379.  
  1380.     The RSA public key cryptosystem is patented by the Massachusetts
  1381.     Institute of Technology (U.S. patent #4,405,829).  This patent does
  1382.     not apply outside the USA.  Public Key Partners (PKP) holds the
  1383.     exclusive commercial license to sell and sub-license the RSA public
  1384.     key cryptosystem.  The author of this software implementation of the
  1385.     RSA algorithm is providing this software for educational use only. 
  1386.     Licensing this algorithm from PKP is the responsibility of you, the
  1387.     user, not Philip Zimmermann, the author of this software.  The author
  1388.     assumes no liability for any breach of patent law resulting from the
  1389.     unlicensed use of this software by the user.
  1390. *********************************************************************/
  1391.  
  1392.  
  1393. int rsa_decrypt(unitptr M, unitptr C,
  1394.     unitptr d, unitptr p, unitptr q, unitptr u)
  1395.     /*    Uses Quisquater's Chinese Remainder Theorem shortcut 
  1396.         for RSA decryption.
  1397.         M is the output plaintext message.
  1398.         C is the input ciphertext message.
  1399.         d is the secret decryption exponent.
  1400.         p and q are the secret prime factors of n.
  1401.         u is the multiplicative inverse of p, mod q.
  1402.         Note that u is precomputed on the assumption that p<q.
  1403.         n, the common modulus, is not used here because of the
  1404.         Chinese Remainder Theorem shortcut.
  1405.     */
  1406. {
  1407.     unit p2[MAX_UNIT_PRECISION];
  1408.     unit q2[MAX_UNIT_PRECISION];
  1409.     unit temp1[MAX_UNIT_PRECISION];
  1410.     unit temp2[MAX_UNIT_PRECISION];
  1411.     int status;
  1412.  
  1413.     mp_init(M,1);    /* initialize result in case of error */
  1414.  
  1415.     if (mp_compare(p,q) >= 0)    /* ensure that p<q */
  1416.     {    /* swap the pointers p and q */
  1417.         unitptr t;
  1418.         t = p;  p = q; q = t;
  1419.     }
  1420.  
  1421. /*    Rather than decrypting by computing modexp with full mod n
  1422.     precision, compute a shorter modexp with mod p precision... */
  1423.  
  1424. /*    p2 = [ (C mod p)**( d mod (p-1) ) ] mod p        */
  1425.     mp_move(temp1,p);
  1426.     mp_dec(temp1);            /* temp1 = p-1 */
  1427.     mp_mod(temp2,d,temp1);    /* temp2 = d mod (p-1) ) */
  1428.     mp_mod(temp1,C,p);        /* temp1 = C mod p  */
  1429.     status = mp_modexp(p2,temp1,temp2,p);
  1430.     if (status < 0)    /* mp_modexp returned an error. */
  1431.         return(status);    /* error return */
  1432.  
  1433. /*    Now compute a short modexp with mod q precision... */
  1434.  
  1435. /*    q2 = [ (C mod q)**( d mod (q-1) ) ] mod q        */
  1436.     mp_move(temp1,q);
  1437.     mp_dec(temp1);            /* temp1 = q-1 */
  1438.     mp_mod(temp2,d,temp1);    /* temp2 = d mod (q-1) */
  1439.     mp_mod(temp1,C,q);        /* temp1 = C mod q  */
  1440.     status = mp_modexp(q2,temp1,temp2,q);
  1441.     if (status < 0)    /* mp_modexp returned an error. */
  1442.         return(status);    /* error return */
  1443.  
  1444. /*    Now use the multiplicative inverse u to glue together the
  1445.     two halves, saving a lot of time by avoiding a full mod n
  1446.     exponentiation.  At key generation time, u was computed
  1447.     with the secret key as the multiplicative inverse of
  1448.     p, mod q such that:  (p*u) mod q = 1, assuming p<q.
  1449. */
  1450.     if (mp_compare(p2,q2) == 0)    /* only happens if C<p */
  1451.         mp_move(M,p2);
  1452.     else
  1453.     {    /*    q2 = q2 - p2;  if q2<0 then q2 = q2 + q  */
  1454.         if (mp_sub(q2,p2))    /* if the result went negative... */
  1455.             mp_add(q2,q);        /* add q to q2 */
  1456.  
  1457.         /*    M = p2 + ( p * [(q2*u) mod q] )     */
  1458.         mp_mult(temp1,q2,u);        /* temp1 = q2*u  */
  1459.         mp_mod(temp2,temp1,q);        /* temp2 = temp1 mod q */
  1460.         mp_mult(temp1,p,temp2);    /* temp1 = p * temp2 */
  1461.         mp_add(temp1,p2);        /* temp1 = temp1 + p2 */
  1462.         mp_move(M,temp1);        /* M = temp1 */
  1463.     }
  1464.  
  1465.     mp_burn(p2);    /* burn the evidence on the stack...*/
  1466.     mp_burn(q2);
  1467.     mp_burn(temp1);
  1468.     mp_burn(temp2);
  1469.     return(0);    /* normal return */
  1470. }    /* rsa_decrypt */
  1471.  
  1472.  
  1473. /****************** end of MPI library ****************************/
  1474.